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Abstract 

b : 

The results of 3D modelling of the flow structure in the classical 
symbiotic system Z Andromedae are presented. Outbursts in systems 
of this type occur when the accretion rate exceeds the upper limit of 
the steady burning range. Therefore, in order to realize the transition 
from a quiescent to an active state it is necessary to find a mechanism 
able to sufficiently increase the accretion rate on a time scale typical 
to the duration of outburst development. 

Our calculations have confirmed the transition mechanism from 
quiescence to outburst in classic symbiotic systems suggested earlier 
on the basis of 2D calculations (Bisikalo et al, 2002). The analysis of 
our results have shown that for wind velocity of 20 km/s an accretion 
disc forms in the system. The accretion rate for the solution with 
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the disc is ~ 22.5 — 25% of the mass loss rate of the donor, that 
is, ~ 4.5 — 5 • 10" 8 M Q /yr for Z And. This value is in agreement with 
the steady burning range for white dwarf masses typically accepted 
for this system. When the wind velocity increases from 20 to 30 km/s 
the accretion disc is destroyed and the matter of the disc falls onto 
the accretor's surface. This process is followed by an approximately 
twofold accretion rate jump. The resulting accretion rate growth is 
sufficient for passing the upper limit of the steady burning range, 
thereby bringing the system into an active state. The time during 
which the accretion rate is above the steady burning value is in a very 
good agreement with observations. 

The analysis of the results presented here allows us to conclude that 
small variations in the donor's wind velocity can lead to the transition 
from the disc accretion to the wind accretion and, as a consequence, 
to the transition from quiescent to active state in classic symbiotic 
stars. 

1 Introduction 

Symbiotic stars are characterized by peculiar spectra where molecular absorp- 
tion bands - the characteristic features of a cool giant - are present together 
with emission lines corresponding to a high excitation level. The red and 
infrared spectra of symbiotic stars are typical of cool giants while in the UV 
range they are characterized by a very hot continuum. It is widely assumed 
that most symbiotic stars are detached binaries consisting of a cool giant 
and a white dwarf surrounded by a nebulosity 0. Mass exchange in these 
systems is driven by stellar wind. 

The goal of this work is to carry out a study of classical symbiotic systems 
using as example one of the most well-studied representatives of this class - 
Z And. According to the analyses of its energy distribution in a wide spectral 
range j2], the components of Z And have the following characteristics: the 
cool M3.5III giant has a mass ~ 2 M & and a radius ~ 100-R©, while the 
white dwarf has a mass ~ 0.6 M & , a radius ~ O.O7-R and a temperature 
~ 10 5 K. The giant loses mass at the rate of ~ 2 x 10~ 7 M Q /yr. Gas in 
the circumbinary envelope has an electron density of ~ 2 • 10 10 cm" 3 and a 
temperature ~ 1.5 -v- 8 ■ 10 4 K. The separation of the system is 482i? Q , the 
orbital period is 758 days. More than 100 years of observations have shown 
that this star presents outbursts. The last of them took place in 2000 and 
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has been actively observed at different wavelengths. 

Following the characteristics of their outbursts, symbiotic stars can be di- 
vided in two types [3J 0] . The first type includes so-called classical symbiotics 
(e.g. Z And, AG Peg, AG Dra, CI Cyg, AX Per). Their outbursts typically 
last a few months and have amplitudes of 2-4 m . The typical feature of such 
outbursts is the weakening of high excitation lines while the brightness in- 
creases. The detailed description of the behaviour of these stars is given in 
1]. Symbiotic novae (e.g. V1016 Cyg, V1329 Cyg, RS Oph, HM Sge) show 
the second type of outbursts. These outbursts are longer and brighter. Their 
energetics can even exceed that of novae outbursts. The principal difference 
of these outbursts from those of the first type is the increase of the ionization 
degree while the brightness increases. The behaviour of these stars during 
outbursts is described in |H1 19110] . The second type of outbursts is similar 
to slow classical novae outbursts [T""] D3 (for differences between symbiotic 
novae and slow classical novae outbursts see "13]). 

The most probable mechanism describing the observational manifesta- 
tions of classical symbiotics as well as of symbiotic novae is thermonuclear 
burning on the accretor's surface [see, e.g. K*] [13] . It has been shown in var- 
ious works [TU HS1 EJ [T"""[ E] that the process of thermonuclear burning of 
hydrogen on the white dwarf surface strongly depends on the accretion rate 

M accr [3J. Only, in a narrow range of M accr is the steady hydrogen burning 
possible "T*"l [TBI fi"*"] . The lower limit of this range is given by the expression 

prong 

/ M \ 3 ' 57 

M steady ,min = 1-3 • 10~ 7 • ( — J M /yr , 

while the higher limit is given by [22J 12*3] 

M steady , max = 6 ■ 1(T 7 ■ (J£- - 0.522^) M /yr , 

For Z And where M = O.6M , steady burning is possible in the range 2.1 • 

l(T 8 M /yr < M accr < 4.7 • lO~ 8 M /yr. If the accretion rate is below this 
interval, thermonuclear burning exhausts the matter faster than accretion 
can replenish it bringing the burning to a halt. Following this, hydrogen 
begins to accumulate on the white dwarf surface until the pressure near the 
hydrogen envelope base reaches a critical value and a hydrogen shell flash 
occurs. The bolometric luminosity increases by a factor 10-100 times in the 
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course of ~ 1 year and stays at the active state for ~ 10 years. The total 
duration of the outburst is dozens of years that corresponds to the nuclear 
time scale. This is observed in symbiotic novae. 

If for some reason M accr becomes larger than M steady, max, the accreted 
matter accumulates above the burning shell and expands to giant dimensions. 
Such behaviour is typical for classical symbiotic stars. Though the bolometric 
luminosity stays constant, the visual brightness increases by 1-2 magnitudes 
and the effective temperature decreases. Optical outburst develops on the 
thermal timescale. In accordance with existing models of thermonuclear 
burning on a white dwarf surface |2()|15| . small variations of the accretion 
rate can lead to significant changes of the temperature at practically constant 
bolometric luminosity. For a white dwarf of 1.2M mass, an envelope mass 
change of 4 x 10~ 7 M & can result in a 3 m luminosity increase [3]. This value 
corresponds to the typical amplitude of Z And outbursts. 

Z And belongs to the class of classical symbiotics and, in accordance with 
the accepted model, the accretion rate in the quiescent state should be in the 
steady burning range. For an outburst to develop the accretion rate should 
become large enough to reach a value above the steady burning range. It 
should be noted that the time scale of the accretion rate increase should 
correspond to the characteristic time of outburst development (~ 100 days). 
Such behaviour of the accretion rate in a binary system is possible in the 
framework of the mechanism proposed in using results of 2D gasdynamic 
modelling. The main idea of the mechanism is that even minor variations in 
the donor's wind velocity is enough to lead to the accretion regime change - 
from disc accretion to wind accretion. During the transition period, namely 
during the disc destruction the accretion rate increases abruptly and exceeds 
the upper limit of the steady burning range. The goal of this work is to 
check the possibility of realizing such a mechanism using a more realistic 3D 
model. 

2 The model 

In order to study the gas flow structure in the symbiotic system Z And 3D 
numerical simulations have been carried out. The zero point of the coordinate 
system was placed at the centre of the accretor, the x- axis was directed along 
the line connecting centres of the components and oppositely to the mass- 
losing star, the y-axis - along the accretor orbital motion, and the z-axis - 
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along the axis of rotation of the binary system. The flow was described by 
the system of Euler equations in corotating coordinate frame: 

dp dpu dpv dpw 
dt dx dy dz 



dpu d(pu 2 + P) dpuv dpuw <9$ „ _ 
+ ^-^ + -7T- + -T- = -p-z- + 2Qvp. 



dt dx dy dz dx 



dpv ^ dpuv ^ d(pv 2 + -P) _|_ dpvw <9$ ^ 
dt dx dy dz dy 



dpw dpuw dpvw d(pw 2 + P) <9$ 
dt dx dy dz ^ dz 



dpE dpuh dpvh dpwh <9$ <9$ <9$ 

dt dx dy dz ^ dx ^ dy ^ dz 

Here u = (it, v, w) is the velocity vector, P - the pressure, p - the density, 
h = e + P/p + \u\ 2 /2 - the specific total enthalpy, E = e + \u\ 2 /2 - the 
specific total energy, e - the specific internal energy, Q - the angular velocity 
of binary system's rotation, and $(r) - the force potential. 

In the standard definition when only the gravitational forces from the 
point mass components and the centrifugal force are taken into account, the 
force potential is given by: 



_ ( GMi GM 2 x ., , 2 

$(r) = - 1 - - y 2 Jr(r - r c y 

- rd \r - r 2 \ 



Here Mi is the mass of the accretor, M 2 - the mass-losing star's mass, r 1; r 2 - 
the radius-vectors of the centres of components, r c - the radius-vector of the 
centre of mass of the system. This is the so-called Roche potential. But in our 
case the additional force responsible for the donor's wind acceleration should 
also be taken into account. Therefore, the form of the potential changes. 
Previous studies (e.g., have shown that the general flow structure 

in the system where components do not fill their Roche lobes is defined first 
and foremost by the stellar wind parameters. Since the mechanism of gas 
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acceleration is poorly known for cool giants, we mimic the radiation pressure 
by reducing the gravitational attraction force to the donor star by a factor 
(1 — T). So, the modified force potential looks the following: 



$(r) = -(i - r)^ — — - j , - y 2 n 2 (r - r c ) 2 

r 2\ 



To close the system the perfect gas equation of state was used. 

P=( 1 -l)pe, 

The value of adiabatic index has been accepted 7 = 1.01, that corresponds 
to the case close to the isothermal one |2TH2T?] . 

We adopted for the parameters of the binary system those of Z And. In 
our numerical simulations all variables were put in dimensionless form as 
follows: the length was normalized by the separation of two stars, A, and the 
time by Q" 1 , and so the velocity was normalized by AQ. The density was 
normalized by that at the surface of the mass-losing star. 

The method of computational fluid dynamics is described in [HO] and 
details are explained in [3*T| We use the simplified flux vector splitting 
(SFS) finite volume scheme to discretize the Euler equation (the description 
of the SFS scheme is given in the appendix of [HI]). 

We assume a symmetry about the orbital plane, and therefore only com- 
pute the upper half-space of the computational domain, which is —2 < x < 1, 
-1.5 < y < 1.5, < z < 1.0. The region is divided into 253 x 253 x 85 
cells. The inner numerical boundary surrounding the mass-losing star is rep- 
resented by an equipotential surface with a mean radius R 2n (we set the 
numerical boundary at a little larger radius, R 2n = lOli? , rather than R 2 ). 
The radius of the mass-accreting star is too small to fit in our numerical grid 
so it is represented by one cell. 

The cells just inside of the inner numerical boundary about the mass- 
losing star are assumed to be filled with gas with dimensionless density p 2 = 
1, dimensionless sound speed c s = 0.178 (corresponding to T 2 = 3200-fT), and 
a dimensionless normal speed to the surface, V. The real normal velocity 
of the gas on the numerical boundary is determined by solving Riemann 
problems between the cells adjacent to the inner boundary surface. The cell 
representing the mass-accreting star and the cells just outside of the outer 
numerical boundary are filled by a gas with density p = 10~ 9 , pressure 
p = 10~ 8 /7, and three velocity components u = v = w = 0. This 



7 



assumption does not mean that accretion onto the mass accreting star nor 
that escape from the computational domain do not occur. The velocity of 
the gas on the boundaries are computed by solving Riemann problems. At 
the initial stage t = 0, the computational domain except the region inside 
of the inner boundaries is assumed to be filled by a gas with p = 10~ 9 , 
p = 10~ 8 /7, and u = v = w = 0. This gas is gradually replaced by the 
gas supplied from the mass-losing star. We compute up to a few rotation 
periods that is generally long enough for the system to reach a quasi-steady 
state. 



3 Results 

The velocity of the mass- losing star's wind in Z And is ~25 km/s 0. It 
is presumed in the framework of the considered mechanism that for slightly 
smaller velocities an accretion disc will form while for slightly greater values 
the disc will be destroyed and a wind accretion type of flow will occur. In 
order to check the applicability of the proposed mechanism, numerical sim- 
ulations for velocities of 25±5km/s have been carried out. 1 Calculations for 
the wind velocity V w = 20 km/s were necessary for studying the possibility 
of disc formation in the system, accretion rate estimation and to check if the 
accretion rate is indeed within the limits of the steady burning range. Cal- 
culations for V w = 30 km/s were conducted in order to obtain parameters of 
the flow structure in the system without disc. To answer the questions con- 
cerning processes taking place after the wind velocity change - if the disc will 
be destroyed, how the accretion rate will grow, if the accretion rate change 
will go above the upper limit of the steady burning range - calculations with 
a wind velocity jump from 20 — > 30 km/s have been carried out Fig.l 

The results of our numerical modelling have shown that for wind velocity . . 

Fia 2 

V w = 20 km/s and values of the parameter T 6 [0.85; 0.95] a steady accretion I y ' 

disc forms in the system. In Figures 1,2 the flow structure for the case 
V w = 20 km/s and T=0.94 is shown. Density contours and velocity vectors in 
the equatorial plane of the system are presented for the whole computational 
domain (Fig. 1) and for the near-accretor area [— 0.5v4 . . . A] x [—A . . . 0.5 A] 
(Fig. 2). The results of the calculations are presented at the time t ~ 5P or b 
when a steady regime was already established in the system. These results 

Small variations of the wind velocity within ±5 km/s can be easily explained by an 
activity of the giant @] . 
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show that for this wind velocity value a bow shock and an accretion disc of 
approximately 50 — 60R & form in the system. Note, that two spiral density 
waves are seen in the accretion disc. In this solution the accretion rate is 
~ 22.5 — 25% of matter lost by the donor. Similar values of the accretion 
efficiency for solution with accretion disc were obtained by Mastrodemos 
and Morris for detached binaries. For Z And where the mass loss rate is 
estimated to be ~ 2 ■ 10~ 7 M /yr, the accretion rate will be ~ 4.5 — 5 • 10 -8 
Mq/jy, i.e. corresponding to the steady burning range (near the upper limit). 

As one can see, due to the influence of the pressure gradient and the pres- 
ence of spiral shocks, the disc structure is far from a conventional Keplerian 
disc. Since we used the full set of Euler equations that incorporate advec- 
tion term in energy equation, the obtained flow structure is more similar to 
advection-dominated discs as discussed by Walder and Folini pMj . Of course, 
the system of Euler equations with adiabatic energy equation does not in- 
corporate viscous heating nor radiative cooling, but the numerical viscosity, 
on the one hand, and the choice of the adiabatic index 7 close to 1, on the 
other hand, provide an implicit account of these processes. Thus, we believe 
that our mathematical model is adequate to the physics of accretion discs. Fig. 3 

The solution with wind velocity V w = 30 km/s is presented in Figures 3 , — : — 
and 4. Just as in Figures 1,2, here density contours and velocity vectors are I ^9-4 
presented for all the computational domain (Fig. 3) and for the area near 
the accretor [— 0.5 A . . . A] x [—A . . . 0.5A] (Fig. 4). The situation presented 
corresponds to the moment when the steady regime has already become 
established. The analysis of these results show that in the case when wind 
velocity equals 30 km/s the cone shock is close to the accretor and does not 
leave any room for the formation of a disc. In this solution a wind accretion 
type of flow instead of disc accretion takes place with the accretion rate being 
~ 11 — 13% of the matter leaving donor's surface. Similar flow structures were 
obtained in 3D simulations performed by Dumm et al. jHSl- According to 
their results, the accretion efficiency equals 6% for the case of wind accretion. Fig. 5 

The analysis of results presented above allows to conclude that even small 
change in the donor's wind velocity (within ±5 km/s around the observed 
value 25 km/s) leads to a flow structure and accretion regime changeover, 
namely, to the transition from a disc to a wind accretion flow. In order to 
consider the process of transition between these two regimes after the wind 
velocity increase, we took the stationary solution for the case 20 km/s after 
~ 5P or b from the beginning of calculations (situation presented in Figures 
la, lb) and then raised the velocity up to 30 km/s. After the increase of 
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the wind velocity on the inner boundary the flow structure changes, evolving 
from a state with an accretion disc to the one with a cone shock presented in 
Figures 2a, 2b. Let us consider the process of the flow restructuring. After 
0.17P or b (~133 days) since the wind velocity change, the matter moving with 
increased velocity reaches the vicinity of the accretor, namely, the bow shock 
located in front of it. Then the wind continues to move further crushing the 
accretion disc and making the matter of the disc fall on the accretor's surface. 
A snapshot of this flow rearrangement is presented in Figure 5. The area size 
and all the designations are the same as in Figures 2,4. This situation corre- 
sponds to ~180 days after the wind velocity increase. It should be mentioned 
that the study of the transition period is limited in the framework of this 
model. After the accretion rate jump the accepted boundary conditions on 
the accretor change and the model used doesn't describe the real physical 
situation any more. Correspondingly, the presented results of calculations of 
the flow rearrangement period are correct only at first stages. 

The behaviour of the accretion rate with time is shown in Figures 4a,b. 
Here the time passed from the beginning of the calculation in units of t = 
P or b/2ir is shown on the x axis. The accretion rate is given in dimensionless 
units. 2 From the moment when the matter with increased velocity reaches 
the vicinity of the accreting component, the accretion rate begins to grow 
and reaches its maximum after approximately 0.06P or b (~47 days). The 
maximum value of the accretion rate is 2-2.2 times as compared to the state 
with disc accretion. 

As shown above, for accretor of O.6M mass, the ratio of the upper limit 

of the steady burning range to its lower limit M steady , max / M steady , min w 2.2. 
For a white dwarf of O.55M mass (the value given in 0]) an increase of 
10% only should be enough to exceed the upper limit. So, the approximately 
twofold accretion rate growth obtained in our calculations is large enough to 
transfer the system from the quiescent to active state. 

The analysis of the data presented in Figures 4a,b shows that the time of 
full disc destruction is ~ 180 days. The time during which the accretion rate 
exceeds the upper limit of the steady burning range should correspond to 
the time of outburst development. If we suppose that the exceedance of the 
limit of the steady burning interval occurs when the accretion rate increases 
by 1.5 times, the time during which accretion rate is above this value will 

2 In these units the rate of the donor's mass loss is equal to 0.355 for the solution with 
a disc and 0.535 for the solution with a cone shock. 
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be approximately 100 days and ~ 2 x 10 M & will be accreted during this 
time. 



4 Conclusions 

Outbursts in classical symbiotic stars occur when the accretion rate becomes 
greater than the upper limit of the steady burning range. For transition from 
quiescent to active state it is necessary to provide an increase in the accretion 
rate on a rather long time interval corresponding to the characteristic time 
of the outburst development (~100 days). Earlier, in |25> a mechanism pro- 
viding the required accretion rate increase in the system has been proposed 
on the basis of 2D calculations. According to this mechanism even minor 
change of the donor's wind velocity is enough for the accretion regime to 
change. During the transition from disc to wind accretion the accretion disc 
is destroyed and the wind with increased velocity makes the matter of the disc 
fall onto the accretor's surface. The analysis of 2D calculations has shown 
that during this process the accretion rate growth is enough for system to 
overstep the limits of steady burning, therefore, leading to the development 
of an outburst. 

In order to answer if this mechanism can work in observed astrophysi- 
cal objects, the gasdynamic modelling of the flow structure in the classical 
symbiotic system Z And using more realistic 3D model has been carried out. 
The results of our calculations allowed us to draw following conclusions: 

1. It was found that for the donor wind velocity 20 km/s and T = 0.85 — 
0.95 the accretion disc forms in the system. The accretion rate for the 
solution with disc is ~ 22.5 — 25% of the matter lost by the donor. 

2. If the wind velocity equals 30 km/s the accretion disc disappears and 
a cone shock forms. The accretion rate for this case is ~ 11 — 13% of 
the matter that left the donor. 

3. These two solutions show that even minor change in the donor's wind 
velocity (within the limits of ±5km/s from the observed value 25 km/s 
j2j) leads to the change in flow structure and to the accretion regime 
change, namely, to the transition from the disc accretion to the wind 
accretion. 
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4. In accordance with the proposed mechanism, the solution with the 
accretion disc should correspond to a quiescent state of the system. 
According to our results, the accretion rate in this case is ~ 22.5 — 
25%. So, for Z And, where the mass loss rate is estimated to be ~ 
2 ■ lO~ 7 M /yr, the accretion rate will be ~ 4.5 — 5 ■ lO~ 8 M /yr, in 
good agreement with the steady burning interval for the most probable 
white dwarf masses. 

5. In the solution where the wind velocity increases from 20 to 30 km/s 
the accretion disc is destroyed and the matter it contained falls onto 
the accretor's surface. 

6. The detailed study of the accretion regime change after the wind veloc- 
ity increase has shown that the process of disc destruction is followed 
by an accretion rate jump. The value of maximum accretion rate is 
~2-2.2 times as much as the initial value. As it has been mentioned 
above, the steady burning range is rather narrow and if the accretor's 
mass equals O.6M the increase of the accretion rate by 2.2 times is 
unambiguously enough for overstepping the upper limit of this range. 
Results of our calculations give an accretion rate increase large enough 
to put the system outside the steady burning range and to move it up 
into the active state. 

7. The time during which the accretion rate is above the steady burning 
range should correspond to the outburst development time, that is 
approximately 100 days for Z And. According to our results the full 
time of disc destruction is approximately ~ 180 days. If we assume 
that the system leaves the steady burning range after the accretion 
rate increases in 1.5 times, the corresponding time will be ~ 100 days, 
in good agreement with observations. 

8. The typical amplitude of Z And outburst is ~ 3 m . For a white dwarf of 
1.2M mass such brightness increase is provided by an envelope mass 
~ 4 x 1O~ 7 M In our calculations the envelope mass is smaller, 
~ 1O _8 M . 3 It is obvious, that for correct comparison of envelope 
masses the calculations of thermonuclear burning on the surface of a 



3 It should be noted that in the solution obtained the mass of the accreted envelope 
depends on the model parameters (wind velocity, disc temperature and density etc.) and 
is, therefore, not a fixed value. 
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white dwarf with the mass equal to that of the accretor in Z And are 
required. Moreover it is necessary to take into account the asymmetry 
of accretion. Unfortunately, we do not know any works where such 
estimations were made, so, the question on the envelope mass remains 
open. Given the accretion rate corresponding to the upper limit of the 
steady burning range and given the typical time of the outburst of ~ 
1/3 year, the envelope mass should not exceed significantly ~ 1O~ 8 M , 
in good agreement with our results. 

The analysis of results presented here allows us to draw the main con- 
clusion that the transition from quiescent to active state in symbiotic stars 
can be concerned with the accretion regime change (transition from the disc 
accretion to wind accretion) as the result of insignificant variations of donor's 
wind velocity. 
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FIGURE CAPTIONS 
for the paper by Bisikalo et al "3D Gasdynamic Modelling ..." 

Fig. 1. Density contours and velocity vectors in the equatorial plane of the 
system for the case V w = 20 km/s and T=0.94. The empty circle centered at 
(-1,0) corresponds to the donor (radius of the circle equals donor's radius), 
the diamond at (0,0) point marks the accretor. 

Fig. 2. The same as in the Fig. \I\m the vicinity of the accretor. 

Fig. 3. Density contours and velocity vectors in the equatorial plane of the 
system for the case V w = 30 km/s and T=0.94. All designations are the same 
as in Fig. [T] 

Fig. 4. The same as in Fig. El in the vicinity of the accretor. 

Fig. 5. Density contours and velocity vectors in the area near the accretor 
~180 days after the wind velocity increase from 20 to 30 km/s. 

Fig. 6. Accretion rate change for the solution where the wind velocity is 
increased from 20 to 30 km/s. The vertical line marks the moment of the 
wind velocity change. 

Fig. 7. The same as in Fig. |H1 for short time interval beginning from the 
moment of wind velocity increase from 20 to 30 km/s. 
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Figure 4: For the paper by Mitsumoto et al. "3D Gasdynamic Modelling . . ." 
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Figure 5: For the paper by Mitsumoto et al. "3D Gasdynamic Modelling . . ." 
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Figure 6: For the paper by Mitsumoto et al. "3D Gasdynamic Modelling ..." 
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Figure 7: For the paper by Mitsumoto et al. "3D Gasdynamic Modelling . . . 



